!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine weno_lx (fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen, fxpo,                                                          &
   ro, mx, my, mz, bx, by, bz, en, po, nx, ny, gm, ch2, riemann_solver_flag)
!==================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx, ny, riemann_solver_flag

   real(8), intent(in) :: gm, ch2
   real(8), dimension(nx, ny), intent(in) :: ro, mx, my, mz, bx, by, bz, en, po
   real(8), dimension(nx, ny), intent(inout) :: fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen, fxpo

   integer(4) :: i, ip1, im1, j

   real(8), dimension(nx, ny) :: ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, po_l
   real(8), dimension(nx, ny) :: ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, po_r
!   real(8), dimension(nx, ny) :: s

   real(8) :: fxro_l, fxmx_l, fxmy_l, fxmz_l, fxbx_l, fxby_l, fxbz_l, fxen_l, fxpo_l
   real(8) :: fxro_r, fxmx_r, fxmy_r, fxmz_r, fxbx_r, fxby_r, fxbz_r, fxen_r, fxpo_r
   real(8) :: alpha_0, alpha_1, d_0, d_1, var_0, var_1
   real(8) :: eps, gmm1, vx, vy, vz, b2, v2, pr, max_speed, c2, s2, ca2, cfx

!----------------------------------------------------------------------------------------------------------------------------------|
   eps = 1.0d-12
   gmm1 = gm - 1.0d0

   d_0 = 2.0d0 / 3.0d0
   d_1 = 1.0d0 / 3.0d0
   do j = 1, ny
   do i = 2, nx - 1
      ip1 = i + 1
      im1 = i - 1
      var_0 = 0.5d0 * (ro(i, j) + ro(ip1, j))
      var_1 = 0.5d0 * (3.0d0 * ro(i, j) - ro(im1, j))
      alpha_0 = d_0 / ((eps + (ro(ip1, j) - ro(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (ro(i, j) - ro(im1, j)) ** 2) ** 2)
      ro_l(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (mx(i, j) + mx(ip1, j))
      var_1 = 0.5d0 * (3.0d0 * mx(i, j) - mx(im1, j))
      alpha_0 = d_0 / ((eps + (mx(ip1, j) - mx(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (mx(i, j) - mx(im1, j)) ** 2) ** 2)
      mx_l(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (my(i, j) + my(ip1, j))
      var_1 = 0.5d0 * (3.0d0 * my(i, j) - my(im1, j))
      alpha_0 = d_0 / ((eps + (my(ip1, j) - my(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (my(i, j) - my(im1, j)) ** 2) ** 2)
      my_l(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (mz(i, j) + mz(ip1, j))
      var_1 = 0.5d0 * (3.0d0 * mz(i, j) - mz(im1, j))
      alpha_0 = d_0 / ((eps + (mz(ip1, j) - mz(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (mz(i, j) - mz(im1, j)) ** 2) ** 2)
      mz_l(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (bx(i, j) + bx(ip1, j))
      var_1 = 0.5d0 * (3.0d0 * bx(i, j) - bx(im1, j))
      alpha_0 = d_0 / ((eps + (bx(ip1, j) - bx(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (bx(i, j) - bx(im1, j)) ** 2) ** 2)
      bx_l(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (by(i, j) + by(ip1, j))
      var_1 = 0.5d0 * (3.0d0 * by(i, j) - by(im1, j))
      alpha_0 = d_0 / ((eps + (by(ip1, j) - by(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (by(i, j) - by(im1, j)) ** 2) ** 2)
      by_l(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (bz(i, j) + bz(ip1, j))
      var_1 = 0.5d0 * (3.0d0 * bz(i, j) - bz(im1, j))
      alpha_0 = d_0 / ((eps + (bz(ip1, j) - bz(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (bz(i, j) - bz(im1, j)) ** 2) ** 2)
      bz_l(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (en(i, j) + en(ip1, j))
      var_1 = 0.5d0 * (3.0d0 * en(i, j) - en(im1, j))
      alpha_0 = d_0 / ((eps + (en(ip1, j) - en(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (en(i, j) - en(im1, j)) ** 2) ** 2)
      en_l(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (po(i, j) + po(ip1, j))
      var_1 = 0.5d0 * (3.0d0 * po(i, j) - po(im1, j))
      alpha_0 = d_0 / ((eps + (po(ip1, j) - po(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (po(i, j) - po(im1, j)) ** 2) ** 2)
      po_l(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)
   enddo
   enddo

   d_0 = 1.0d0 / 3.0d0
   d_1 = 2.0d0 / 3.0d0
   do j = 1, ny
   do i = 2, nx - 1
      ip1 = i + 1
      im1 = i - 1
      var_0 = 0.5d0 * (3.0d0 * ro(i, j) - ro(ip1, j))
      var_1 = 0.5d0 * (ro(i, j) + ro(im1, j))
      alpha_0 = d_0 / ((eps + (ro(ip1, j) - ro(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (ro(i, j) - ro(im1, j)) ** 2) ** 2)
      ro_r(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (3.0d0 * mx(i, j) - mx(ip1, j))
      var_1 = 0.5d0 * (mx(i, j) + mx(im1, j))
      alpha_0 = d_0 / ((eps + (mx(ip1, j) - mx(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (mx(i, j) - mx(im1, j)) ** 2) ** 2)
      mx_r(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (3.0d0 * my(i, j) - my(ip1, j))
      var_1 = 0.5d0 * (my(i, j) + my(im1, j))
      alpha_0 = d_0 / ((eps + (my(ip1, j) - my(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (my(i, j) - my(im1, j)) ** 2) ** 2)
      my_r(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (3.0d0 * mz(i, j) - mz(ip1, j))
      var_1 = 0.5d0 * (mz(i, j) + mz(im1, j))
      alpha_0 = d_0 / ((eps + (mz(ip1, j) - mz(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (mz(i, j) - mz(im1, j)) ** 2) ** 2)
      mz_r(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (3.0d0 * bx(i, j) - bx(ip1, j))
      var_1 = 0.5d0 * (bx(i, j) + bx(im1, j))
      alpha_0 = d_0 / ((eps + (bx(ip1, j) - bx(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (bx(i, j) - bx(im1, j)) ** 2) ** 2)
      bx_r(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (3.0d0 * by(i, j) - by(ip1, j))
      var_1 = 0.5d0 * (by(i, j) + by(im1, j))
      alpha_0 = d_0 / ((eps + (by(ip1, j) - by(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (by(i, j) - by(im1, j)) ** 2) ** 2)
      by_r(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (3.0d0 * bz(i, j) - bz(ip1, j))
      var_1 = 0.5d0 * (bz(i, j) + bz(im1, j))
      alpha_0 = d_0 / ((eps + (bz(ip1, j) - bz(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (bz(i, j) - bz(im1, j)) ** 2) ** 2)
      bz_r(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (3.0d0 * en(i, j) - en(ip1, j))
      var_1 = 0.5d0 * (en(i, j) + en(im1, j))
      alpha_0 = d_0 / ((eps + (en(ip1, j) - en(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (en(i, j) - en(im1, j)) ** 2) ** 2)
      en_r(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

      var_0 = 0.5d0 * (3.0d0 * po(i, j) - po(ip1, j))
      var_1 = 0.5d0 * (po(i, j) + po(im1, j))
      alpha_0 = d_0 / ((eps + (po(ip1, j) - po(i, j)) ** 2) ** 2)
      alpha_1 = d_1 / ((eps + (po(i, j) - po(im1, j)) ** 2) ** 2)
      po_r(i, j) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)
   enddo
   enddo

   if (riemann_solver_flag .eq. 1) then
      call hlld_x (fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen,                                                                 &
         ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, nx, ny, gm)
   elseif (riemann_solver_flag .eq. 2) then
      call hllc_x (fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen,                                                                 &
         ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, nx, ny, gm)
   elseif (riemann_solver_flag .eq. 3) then
      call roe_x (fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen,                                                                  &
         ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, nx, ny, gm)
   endif

   do j = 1, ny
   do i = 1, nx - 1
      ip1 = i + 1
!      im1 = max (i - 1, 1)
!      max_speed = maxval(s(im1 : ip1, j))

      vx = mx_l(i, j) / ro_l(i, j) 
      vy = my_l(i, j) / ro_l(i, j) 
      vz = mz_l(i, j) / ro_l(i, j) 
      b2 = bx_l(i, j) * bx_l(i, j) + by_l(i, j) * by_l(i, j) + bz_l(i, j) * bz_l(i, j)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_l(i, j) - v2 * ro_l(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1
      c2 = gm * pr
      s2 = c2 + b2
      ca2 = bx_l(i, j) * bx_l(i, j)
      cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_l(i, j) / 2.0d0)
      max_speed = cfx + abs (vx)

      if (riemann_solver_flag .eq. 0) then
         fxro_l = mx_l(i, j)
         fxmx_l = vx * mx_l(i, j) + pr + b2 / 2.0d0 - bx_l(i, j) * bx_l(i, j)
         fxmy_l = vx * my_l(i, j) - bx_l(i, j) * by_l(i, j)
         fxmz_l = vx * mz_l(i, j) - bx_l(i, j) * bz_l(i, j)
         fxbx_l = po_l(i, j)
         fxby_l = vx * by_l(i, j) - bx_l(i, j) * vy
         fxbz_l = vx * bz_l(i, j) - bx_l(i, j) * vz
         fxen_l = (en_l(i, j) + pr + b2 / 2.0d0) * vx - (bx_l(i, j) * vx + by_l(i, j) * vy + bz_l(i, j) * vz) * bx_l(i, j)
      endif
      fxpo_l = ch2 * bx_l(i, j)

      vx = mx_r(ip1, j) / ro_r(ip1, j) 
      vy = my_r(ip1, j) / ro_r(ip1, j) 
      vz = mz_r(ip1, j) / ro_r(ip1, j) 
      b2 = bx_r(ip1, j) * bx_r(ip1, j) + by_r(ip1, j) * by_r(ip1, j) + bz_r(ip1, j) * bz_r(ip1, j)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_r(ip1, j) - v2 * ro_r(ip1, j) / 2.0d0 - b2 / 2.0d0) * gmm1
      c2 = gm * pr
      s2 = c2 + b2
      ca2 = bx_r(ip1, j) * bx_r(ip1, j)
      cfx = sqrt ((s2 + sqrt (s2 * s2 - 4.0d0 * c2 * ca2)) / ro_r(ip1, j) / 2.0d0)
      max_speed = max (cfx + abs (vx), max_speed)

      if (riemann_solver_flag .eq. 0) then
         fxro_r = mx_r(ip1, j)
         fxmx_r = vx * mx_r(ip1, j) + pr + b2 / 2.0d0 - bx_r(ip1, j) * bx_r(ip1, j)
         fxmy_r = vx * my_r(ip1, j) - bx_r(ip1, j) * by_r(ip1, j)
         fxmz_r = vx * mz_r(ip1, j) - bx_r(ip1, j) * bz_r(ip1, j)
         fxbx_r = po_r(ip1, j)
         fxby_r = vx * by_r(ip1, j) - bx_r(ip1, j) * vy
         fxbz_r = vx * bz_r(ip1, j) - bx_r(ip1, j) * vz
         fxen_r = (en_r(ip1, j) + pr + b2 / 2.0d0) * vx - (bx_r(ip1, j) * vx + by_r(ip1, j) * vy + bz_r(ip1, j) * vz) * bx_r(ip1, j)
      endif
      fxpo_r = ch2 * bx_r(ip1, j)

      if (riemann_solver_flag .eq. 0) then
         fxro(i, j) = 0.5d0 * (fxro_l + fxro_r - max_speed * (ro_r(ip1, j) - ro_l(i, j)))
         fxmx(i, j) = 0.5d0 * (fxmx_l + fxmx_r - max_speed * (mx_r(ip1, j) - mx_l(i, j)))
         fxmy(i, j) = 0.5d0 * (fxmy_l + fxmy_r - max_speed * (my_r(ip1, j) - my_l(i, j)))
         fxmz(i, j) = 0.5d0 * (fxmz_l + fxmz_r - max_speed * (mz_r(ip1, j) - mz_l(i, j)))
         fxbx(i, j) = 0.5d0 * (fxbx_l + fxbx_r - max_speed * (bx_r(ip1, j) - bx_l(i, j)))
         fxby(i, j) = 0.5d0 * (fxby_l + fxby_r - max_speed * (by_r(ip1, j) - by_l(i, j)))
         fxbz(i, j) = 0.5d0 * (fxbz_l + fxbz_r - max_speed * (bz_r(ip1, j) - bz_l(i, j)))
         fxen(i, j) = 0.5d0 * (fxen_l + fxen_r - max_speed * (en_r(ip1, j) - en_l(i, j)))
      endif
      fxpo(i, j) = 0.5d0 * (fxpo_l + fxpo_r - max_speed * (po_r(ip1, j) - po_l(i, j)))
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine weno_lx
